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The largest contributions to the n = 2 Lamb-shift, fine structure iaterval and 2s hyperfine structure 
of muonic hydrogen are calculated by exact numerical evaluations of the Dirac equation, rather than 
by a perturbation expansion in powers of 1/c, in the framework of non-relativistic quantum electro- 
dynamics. Previous calculations and the validity of the perturbation expansion for light elements are 
confirmed. The dependence of the various effects on the nuclear size and model are studied. 
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Despite many years of study, the proton charge radius 
has remained relatively poorly known. It has been de- 
rived from measurements in electron-proton collisions 
[1, 2] or from high-precision spectroscopy of hydrogen 
[3-10] as described in the CODATA report in [11]. Tests 
of fundamental physics based on the progress in accu- 
racy of spectroscopy of hydrogen and deuterium have 
been limited by the lack of an accurate value for the pro- 
ton radius. Moreover, the values for the proton radius 
obtained by different methods or different analyses of ex- 
isting experiments are spread over a range larger than the 
uncertainty quoted for the individual results. Two recent 
measurements have resulted in a puzzle. The accurate 
determination of the 2S Lamb shift by laser spectroscopy 
in muonic hydrogen provides a proton size with a ten 
times smaller uncertainty than any previous value and 
it differs hy five standards deviations from the 2006 CO- 
DATA value[12]. At the same time, a new, improved de- 
termination of the charge radius by electron scattering, 
performed at Mainz with the MAMI microtron, provides 
a value in good agreement with the value from hydro- 
gen and deuterium spectroscopy [13, 14]. Taking into 
account improved theory in hydrogen and deuterium 
and the MAMI measurement lead the recently released 
2010 adjustment [15] to differ by 6.9 standard deviation 
between from the proton radius obtained from muonic 
hydrogen. 

Many papers have been published in the last year, 
trying to solve this piizzle. A few are dealing with the 
calculation of the n = 2 level energies in muonic hydro- 
gen. Several others are concerned with the effect of the 
internal structure of the proton on these energies [16-27]. 
Others look at exotic phenomenon beyond the standard 
model [28-34]. 

Many contributions to the Lamb shift, fine, and h5^er- 
fine structure of muonic hydrogen have been evaluated 
over the years; the results are summarized in [35-^0] 
and in a recent book by Eides et al. [41]. Most of these 
calculations are done in the framework of nonrelativis- 
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tic QED. The wavefunction and operators are expanded 
in powers of the fine-structure constant, and the con- 
tributions are obtained by perturbation theory. Hylton 
[42] showed that the perturbation calculation of the fi- 
nite size correction to the vacuum polarization in heavy 
elements gives incorrect results. Since a bound muon 
is closer to the nucleus than a botind electron by a fac- 
tor OT|^/OTe ~ 207, its Bohr radius is slightly smaller than 
the Compton wavelength of the electron Ac - h/rrieC by a 
factor wie/amn ~ 137/207 (a ~ 1/137.036 is the fine struc- 
ture constant, wig and the elecfron and muon mass 
respectively). The Compton wavelength is the scale of 
QED corrections, and for the 2S level, the muon wave- 
function mean radius is only 2.6 times larger than the 
electron Compton wavelength. 

It is thus worthwhile to reconsider the largest correc- 
tions that contribute to the 2S Lamb-shift in muonic hy- 
drogen using non-perturbative methods. In the present 
work, we use the latest version of the MCDF code of 
Desclaux and Indelicato [43], which is designed to cal- 
ciilate properties of exotic atoms [44], to evaluate the 
exact confribution of the elecfron Uehling potential with 
Dirac wavefunctions including the finite nuclear size. In 
the same way, we calculate the Kallen and Sabry contri- 
bution. 

Throughout this paper we will use QED units, H = 1, 
c = 1. The electric charge is given by = Ana. 



I. NUMERICAL EVALUATION OF THE DIRAC 
EQUATION WITH REALISTIC NUCLEAR CHARGE 
DISTRIBUTION MODELS 

A. Evaluation by the numerical solution of the Dirac 
equation 

We calculate higher-order finite size correction, start- 
ing from the Dirac equation with reduced mass, as tech- 
niques for the acciirate niimerical solution of the Dfrac 
equation in a Coulomb potential have been developed 
over a period of many years within theframework of the 
Multiconfiguration Dirac-Fock (MCDF) method for the 
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atomic many-body problem [45-48]. 
The Dirac equation is written as 

[a ■ p + |3fir + Vj^{r)] 0„Kf,(r-) = Gn^.^nK^Xr), (1) 

where a and j3 are the Dirac 4x4 matrices, ^^(r) is the 
Coulomb potential of the nucleus, £„Kfi is the atom total 
energy, and <1> is a one-electron Dirac four-component 
spinor: 



(2) 



in which Xx^tiS^ (p) two-component Pauli spherical 

spinor [46], n is the principal quantum number, k is the 
Dirac quantum number, and /i is the eigenvalue of }z- 
This reduces, for a spherically symmetric potential, to 
the differential equation: 



VNir) ■ 



d_ , K 
dr r 





-£° 


PnAr) ' 


[ QnAr) . 




. QnK(r) . 



(3) 



where PnK{T) and Qn^T) are the large and small radial 
components of the wavefunction, respectively, k the 
Dirac quantum number, E]^^^ is the binding energy, jj.^ 
is the muon reduced mass, f(r = m^Mp/(fn|x -I- Mp) 
and Mp are the muon and proton masses). 

To solve this equation numerically, we use a 5 point 
predictor-corrector method (order h^) [48, 49] on a linear 
mesh defined as 



+ ar„. 



(4) 



with t„ = to + nh, and ro > is the first point of the mesh, 
corresponding to n = 0. This immediately gives to = aro. 
Equation (4) can be inverted to yield 

W(flroe'») 
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where W is the Lambert (or product logarithm) function. 
The wavefunction and differential equation between 
and ro are represented by a 10 term series expansion. For 
a point nucleus, the first point is usually given by ro = 
10~^/Z and h = 0.025. Here we use values down to ro - 
10~^/Z and h = 0.002 to obtain the best possible accuracy. 
For a finite charge distribution, the nuclear boundary is 
fixed at the value r^, where N is large enough to obtain 
sufficient accuracy. The mean value of an operator O, 
that gives the first-order contributions to the energy, is 
calculated as 



Xoo 
dr[Pirf + Q{rf]Oir) 

= JJ dr[P{rf + Q{rf]0{r) 

dt^^[P{rf + Q{rf]0{r) 



(6) 



using 8 and 14 points integration formulas due to 
Roothan. The two integration formulas provide the same 
result within 9 decimal places. 



6. Charge distribution models 

For the proton charge distribution, two models are ex- 
tensively used. The first corresponds to a proton dipole 
(charge) form factor, the second is a gaussian model. 
Here we also use uniform and Fermi charge distribu- 
tiond and fits to experimental data [50, 51]. The analytic 
distributions are parametrized so they provide the same 
mean square radius R. Moments of the charge distribu- 
tion are defined by 



<r">=47T I r^-'"p(r)dr, 
Jo 



(7) 



where the nuclear charge distribution p(r) = pf^{r)/{Ze) 
is normalized by 



p{r)dr = 471 J" p{r)r^dr = 1, 



(8) 



for a spherically S5anmetric charge distribution. The 
mean square radius is R = V< r^ >. 

The potential can be deduced from the charge density 
using tiie well known expression: 

Am r r°° 

Vni^") = I duu^p^{u)-Ane I dMMpN(M). (9) 

Jo Jr 

The exponential charge distribution and correspond- 
ing potential energy are written 



pN(r) = Ze 



Vn{r) = -Ze 



e I 

l-e~i e" 



2a' 



< r" > 



{n + 2)! ET 



(10) 



which gives E, = The gaussian charge distribution 
and potential are given by 



PN(r) = Ze 



V^{r) = -Ze 



e-(f) 

7x3/2^3' 

erf (I) 
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(11) 



where erf is the error function, and E, = -^R. Other 

similar expressions for the two above models can be 
found in [52]. 
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FIG. 1. Electromagnetic interaction between a lepton (narrow 
line) and a nucleon (bold line). 



The electric form factor is related to the charge distri- 
bution by 



Ge(Q') 



dre 



For the exponential model this leads to 
1 



(1 - ¥) 

while for the Gaussian model one has 



(13) 
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(14) 



The two models have an identical slope /6 as functions 
of for q — > as expected (see, e.g, [53]). 

In 1956, Zemach introduced an electromagnetic form 
factor, useful for evaluating the hjrperfine structure en- 
ergy correction 



pem{r) = J p{r- u)ii{u)du, 



(15) 



where ii{u) is the magnetic moment density. Both /.i(m) 
and pern are normalized to unity as in Eq. (8). The 
Zemach radius is given by 



Rz = (rz) 



J rpe„,{r)dr. 



(16) 



The Zemach's radius can be written in momentum space 
as [54, 55] 



-4 r 1 , Gm (Q^) 
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where Kp is the proton anomalous magnetic moment, 
and Gm is normalized so that Gm (0) = 1 + Kp. The expo- 
nential and gaussian models enables to obtain analytic 



results for Rz as a function of the charge and magnetic 
moment radii R and Km- Using (13) or (14) for the expo- 
nential or gaussian model, and Eq. (17), we get respec- 
tively 
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An other useful quantity, which appears in the estima- 
tion of the finite size correction to vacuum polarization 
is the third Zemach's moment 

('''>(2) = / ''V(2)(^)d^, (20) 

where the convolved charge distribution is 



(12) p{2){r) = J p{r - u)p{u)du. 



(21) 



This can be rewritten in the more convenient form [17, 
35], in the limit of large proton masses. 



- 1 + irR' 
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(22) 



It can be easily seen from Eqs. (13) or (14) that the ex- 
pression is finite for q —> 0. 

We now turn to more realistic models, based on ex- 
periment. A recent analysis of the world's data on elas- 
tic electron-proton scattering and calculations of two- 
photon exchange effects provides an analytic expression 
for the electric form factors [51], given as 



Ge(q^) = 1^#^, 



(23) 



where t = q^KAMp). The a, coefficients can be found in 
Table I of Ref. [51]. A second work [50] uses a combi- 
nation of several spectral functions tacking into account 
several resonances and continua like the In, XXand pn 
continua. Here we use the fit resulting from the super- 
convergence approach from this work. This corresponds 
to a sum of 12 dipole-like functions, which are able to 
represent the experimental data with a reduced of 
1.8. A comparison of the electric form factor from both 
works is presented on Fig. 3. It is clear that both ex- 
perimental form factors and the dipole approximation 
with an identical radius are very close. We can obtain 
the charge radius and the next correction by performing 
an expansion in q of the experimental form factors. We 
get for Ref. [51] 



Ge(Q^) 



1 



0.8503^ 



-q^ + 



0.8503^ 
43.3909 



(24) 
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and for Ref . [50] 



0.849952 



0.84995^ 
38.793 



(25) 



These expansions are very close to the one for a dipole 
form factors from Eq. (13). 

In order to compare different charge density models, 
we have performed an analytic evaluation of the charge 
densities corresponding to [51], replacing (23) in (12) and 
performing the inverse Fourier transforms, to obtain the 
corresponding charge distribution, depending on the set 
of a coefficients. The corresponding densities are plot- 
ted and compared to Fermi, Gaussian and exponential 
models. Distances are converted from GeV to fm using 
he = 0.1973269631 GeVfm in the density obtained from 
the experiment. The charge distributions are compared 
on Fig. 2. We plotted both p{r) and p{r)r^ to reveal the 
differences at long and medium distances. The experi- 
mental charge density is rather different from all three 
analytic distributions, while, once multiplied by r^, it is 
closer to the exponential distribution. 




Experiment 



Exponential 




FIG. 2. Top: charge densities p(r), bottom: charge densities 
r^p{r) for the experimental fits in Ref. [51], compared to Gaus- 
sian, Fermi and exponential models distributions. All models 
are calculated to have the same R = 0.850 fm RMS radius as 
deduced from the experimental function. 




FIG. 3. Comparison of the electric form factor from Refs. [50, 
51], withadipole model with the sameR = 0.850 fm as deduced 
from the experimental function. 






FIG. 4. Feynman diagrams corresponding to the full vacuum 
polarization contribution, and expansion in Za. Diagram (a) 
corresponds to the Uehling potential [Eqs. (26) and (28)]. Di- 
agram (b) corresponds to the Wichmann and KroU correction. 
The double line represents a bound lepton wavefunction, the 
wavy line a retarded photon propagator. The single line cor- 
respond to a free electron-positron (or muon-antimuon) pair. 
The grey circles correspond to the interaction with the nucleus. 



II. EVALUATION OF MAIN VACUUM POLARIZATION 
AND FINITE SIZE CORRECTION 



The Fe5mman diagram corresponding to the Uehling 
approximation to the vacuum polarization correction is 
presented in Fig. 4 (a). The evaluation of the vacuum po- 
larization can be performed using standard techniques 
of (perturbative) non-relativistic QED (NRQED) as de- 
scribed in [35, 37]. Here we use the analytic results of 
Klarsfeld [56] as described in [57] and numerical solu- 
tion of the Dirac equation from Sec. I. In order to ob- 
tain higher order effects, we solve the Dirac equation in 
a combined potential resulting from the finite nuclear 
charge distribution and of the Uehling potential. The 
logarithmic singularity of the Uehling potential at the 
origin for a point charge cannot be easily incorporated 
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in a numerical Dirac solver. In the case of a finite charge 
distributions, the singularity is milder, but great care 
must be exercised to obtain results accurate enough for 
oui purpose. 

For a point charge, the Uehling potential, which rep- 
resents the leading contribution to the vacuirai polariza- 
tion, is expressed as [56, 58, 59] 



_ 2a(Za^) l i2 \ 



1.2 1 



-ImpTZ 



(26) 



where nie is the electron mass, Ac is the electron Compton 
wavelength and the function xi belongs to a family of 
functions defined by 

The Uehling potential for a spherically symmetric charge 
distribution is expressed as [56] 

la{Za) 1 r°° , , , , 
yn{r) = 5 dr'r'pir') 

^ ^ Jo 

x[x2(;^|r-r'|)-;i:2(^|r + r'|)](28) 
The expression of the potential at the origin is given by 

Vii(O) = ^ dr'r'p{r')xi [jf], (29) 

while it behaves at large distances as [60] 

(30) 



Vu(r) = 



2a(Za) 1 
371 r 
2 



Xi 
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using the moments of the charge distribution (7). The 
energy shift associated with the potential (26) or (28) in 
first order perturbation is calculated as 



AE";P" = {n, I, K, Hr I Vli I n, I, K, jUr) 



(31) 



where |n, /, k, /ij) is a wavefunction solution of Eq. (1), 
which depends on the reduced mass. It was shown 
recently [61] that this method provides the correct in- 
clusion of the vacuum polarization recoil correction at 
the Barker and Glover level [62]. Since we directly use 
relativistic functions, the other corrections described in 
[61] are automatically included. 



III. HIGHER ORDER QED CORRECTIONS 

A. Reevaluation of the Kallen and Sabry potential 

The Kallen and Sabry potential [63], is a fourth order 
potential, corresponding to the diagrams in Fig. 5. The 






FIG. 5. Feynman diagrams included in the Kallen and Sabry 
V2i(r) potential (Eq. (32)). See Fig. 4 for explanation of sym- 
bols. 



expression for this potential has also been derived on 
Ref. [64-67]. In the previous version of the mdfgme 
code, the Kallen and Sabry potential used was the one 
provided by Ref. [60], which is only accurate to 3 digits. 
The expression of this potential is for a point charge 



(32) 



where 
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(33) 



and 



dx 



(3x2-l)ln(Vx2-l+x) ln(8x(x2-l))^ 



Vx2-1 



(34) 



The function f{t) can be calculated analytically in term 
of the In and dilogarithm functions. Blomqvist [68] has 
shown that Li(r) can be expressed as 



Li{r) = g2{r)]n^r) + gi{r)log{r) + go{r), 



(35) 



and provided a series expansion of this function for small 
r. Fullerton and Rinker [60] provided polynomial ap- 
proximations to the functions gi{r). Here we have nu- 
merically evaluated the function Li(r) to a very good 
accuracy, using Mathematica. We then fitted the coeffi- 
cients of polynomials for the function gi{r). The results 
are presented in Appendix A. For x > 3, we have used 
the functional form 

e~'' (a + b yfr + cr + dr^^^ + ef- -H fr'^^\ 
Lrir) = ^ ^2 -■ (36) 
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The coefficients are also given in Appendix A. 

To obtain the finite nuclear size correction, we use the 
known expression for a spherically-symmetric charge 
distribution [60] 



Vnir) = ^ 



i^\r-r'\) 



^ dr'r'p{r')iLo{ 

-Lo(|-|r + r'|)), (37) 



where 



duLi{u). 



(38) 



Using our approximation to Li{x) in Eqs (35) and (36), we 
obtain the following approximate expressions for Lq{x). 
For X < 3, the expression is very similar to the one for 
Li{x). One obtains 

Loir) = rh2ir) hx^ir) + rh{r)logir) + /io(r), (39) 

The expression for the functions hi are given in Appendix 
B. For the asymptotic function, given for x > 3, we inte- 
grate directly Eq. (36), which yield (fixing the integration 
constant so that Lo(r) is at infinity) 

Loir) = 41.1352787432251923 

- 5.1094977559522696|8.05074798111erf ( ^/r) 



- 1.028091975364Ei(-r) 



+ 



^(-2.02809197536r3/2 + 4.54214815071?- 

■5/2 y 

- 0.494718704003r + 0.98439728916 ^/r 

- 0.344009752879)1 



(40) 



where Ei(r) is the exponential integral. 



IV. NUMERICAL RESULTS 

A. Finite size correction to the Coulomb contribution 

Obtaining the accuracy required from the calculation 
on which has a value of ~ 632.1 eV, while the Lamb 
shift is siO.22 eV with an aim at better than 0.001 meV, 
is a very demanding task. For a point nucleus, we get 
exact degeneracy for the 2s and 2pi/2 Dirac energies. The 
best numerical accuracy was obtained generating the 
wavefunction on a grid with ro = 2 X 10"^ and /i = 2 X 
10"''. This corresponds to ~ 8700 tabulation points for 
the wavefunction, with around 2800 points inside the 
proton. I checked that variations in ro and h do not 
change the final value. The main finite nuclear size effect 
on the 2pi/2-2s energy separation comes from the sum 



of the Dirac energy splitting (the 2pi/2 and 2s level are 
exactly degenerate for a point nucleus). 

I evaluated the different quantities on a grid of proton 
sizes ranging from 0.3 to 1.2 fm, with steps of 0.025 fm (80 
points). I also evaluated the contribution for the muonic 
hydrogen proton size and the COD ATA 2010 proton size. 
The first few terms of the dependence of the relativistic 
energy on the moments of the charge distribution where 
given by Friar [52]. For a Gaussian charge distribution 
he finds for a s state 



(41) 



+ d<r^)<logr) + e<r^)^<logr). 

I use this as a guide to fit my numerical results. 
As a first example a 3-parameter fits provides 

^^2p°"2-2si/2(^^ = -5.19972 + 0.0351289 R' 

-0.0000534235 R*meV (42) 
A better fit is provided by 

AE^^l-_2s,p{R) = -5.19990K^ + 0.03559052?^ 

- 0.000488059R^ + 0.000172334R5 (43) 
-0.0000245051R^meV 

Using Friar functional form with only one logR term, I 
obtain 



ae: 



Coul. 

2pi/2-2si/2 



(R) = -5.199365R2 + 0.03466100R3 
+ 0.00007366037R* 
- 0.00001720960R5 
+ 1.198332 X 10"*^^^ 
+ 0.0002677236R2 log R meV. 



(44) 



The function with two logarithmic terms and close val- 
ues of the BIC and criteria is given by 



AE^;;-2s,,.{R) = -5.199337R2 

+ 0.03458139R^ + 0.0001092856R* 
+ 0.0002788380R2logR 
-0.00004957598R*logR 



(45) 



Criteria for the quality of the fit are plotted in Fig. 6. I 
use both the reduced and a Bayesian information cri- 
terion (BIC) to evaluate the improvement in the value 
when increasing the number of parameters [69]. We ob- 
tain a coefficient for R^ which is -5.1999 meV/fm^ and a 
coefficient for R^ equal to 0.03559 meV/fm^. Using our 
numerical solutions we also find -5.19972 meY/firr and 
0.032908 meV/fm' for the Gaussian model. Borie [39] 
finds -5.1975 meV/fm^ and 0.0347 meV/fm^ for an ex- 
ponential model, and 0.0317 meV/fm^ for a Gaussian 
model, in reasonable agreement with the result pre- 
sented here. 
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shift in Ref. [71]. In the Gaussian model, I find 



FIG. 6. BIC criterium (top) and reduced (bottom) as func- 
tion of the degree of the polynomial and of the logarithmic 
dependence used in the fit of the energy. 



In the same way, I evaluated the finite size correction 
to the fine structure. 



ae: 



Coul. 

2p3/2-2p 



^^^(K) = 8.41563570 
- 0.00005192R2 
+ 1.1818650 xlQ-^K^ 
- 1.19528126 X 10-'*R''meV. 



(46) 



The constant term is in perfect agreement with Borie's 
value 8.41564 meV [70] (Table 7). 

It is interesting to explore at this stage the influence 
of the charge distribution shape on the Coulomb and 
vacuum polarization contribution. Friar and Sick [71] 
have evaluated the third Zemach moment from Eq. (20), 
using the proton-electron scattering data available in 
2005. Using a model-independent analysis, they find 

^r''^^^^ = 2.71 ± 0.13 fm'', leading to an energy shift of 

-0.0247 ± 0.0012 meV. Using the Fourier transform of 
the exponential (13) distribution in Eq. (22), I find 



(2) 



35V3K3 
16 



3.789K^ 



(47) 



showing that proportional to and justifying 

the fit in and P? performed to derive the coefficients 
above. Equation (47) is in exact agreement with the result 
that can be obtained from Eq. (15) in [16], but Eq. (16) 
in the same work is not correct (the denominator should 
be 256, not 64). The value obtained by Friar and Sick 
corresponds to K = 0.894 fm. In that case our energy 
shift is -0.0250 meV in good agreement with the energy 



(2) 



1.960K^ 



(48) 



leading to R = 0.920 fm and an energy shift of 
-0.0256 meV, still in agreement. One can perform a more 
advanced calculation, using the experimental charge 
distribution from Ref. [51], as given in (23). I find 

^r^^ = 2.45 fm^, significantly lower than Friar and Sick's 
value. This lead to R = 0.864 fm for the exponential 
model and R = 0.889 fm for the Gaussian model, provid- 
ing shifts of -0.0226 meV and -0.0232 meV respectively, 
in closer agreement to Borie's work. In a recent paper, De 
Riijula [16] claims that the discrepancy found between 
the charge radii obtained from hydrogen and muonic 
hydrogen could be due to the fact that theoretical calcu- 
lations use too simple a dipole model to represent the 
nucleus. He builds a "toy model" composed of the sum 
of two dipole function corresponding to two resonances 
with different masses. In his model the third moment of 
the charge distribution is much higher than what is de- 
rived from a dipole model, enabling to mostly resolve the 
discrepancy between charge radii obtained from muonic 
and normal hydrogen. He gets 



{r% = 36.6 ± 73 



A3R^ 



(49) 



using R from muonic hydrogen. I use the fit to the exper- 
imental form factor from Ref. [51] as given in Eq.(23) to 
check the result from Ref. [16] against an experimental 
determination. I obtain 



(A ^ 2.4485 ~ 3.98 x 0.850^ 

\ /(2) 



(50) 



very close to the dipole model value of Eq. (47). Using 
the recent MAMI experiment, combined with data from 
[51], Distler et al. [22] obtain 



(r^)^^^ ^ 2.85(8) ~ 4.18 x 0.880^ 



(51) 



The conclusions from Ref. [16], which depend on an 
overly large third moment of the charge distribution are 
thus not supported by experiment. 



B. Finite size correction to the Uehling contribution 

For the vacuum polarization we obtain, for a point 
nucleus. 



AE"'P" , = 205.028201 meV, 



(52) 



to be compared with 205.0282 meV in Ref. [39]. 
Pachucki [37] obtained 205.0243 meV as the sum of the 
non-relativistic 205.0074 meV and first order relativistic 
0.0169 meV corrections. If I calculate the difference be- 
tween (52) and Pachuki non-relativistic value, I obtain a 
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difference of 0.0208076 meV. This is in excellent agree- 
ment with the value provided in Ref. [72], Eq. (6), 
0.020843 meV. 

To achieve this result we used the mesh parameters de- 
scribed in the previous section, and checked by varying 
them so that the results were stable within the dedmal 
places provided here. For finite nuclei, I use the same 
parameters as in the previous section. Again changes in 
ro and h do not change the final value. We get 



AE"'^' , (R) = 205.0282076 - 0.02810970909]?^ 
+ 0.0007111893365K^ 
- 0.00003572368803]?^ 



(53) 



The constant term is in excellent agreement with the one 
in Eq. (52). This result must be combined to Eq. (45) to 
obtain values that can be compared with the literature. I 
obtain: 

AE"*^'*', (R) = 205.0282076 - 5.227446248]?2 
+ 0.03529257801]?^ 
+ 0.00007356191826]?* (54) 
+ 0.0002788380236]?2 log(]?) 
- 0.00004957597920]?* log(]?) 

The R^ coefficient can be compared to the one in Ref. [72] 
Table III, which has the value -5.2254 meVfm"^, which 
contains additional recoil corrections. 

For the Uelhing correction to the fine structiire, I obtain 
in the same way: 

^^J5-2,3/2(^) = 0-0050157881 

- 1.1662334 X IQ-^R^ 
+ 2.2741334 x 10"'']?^ 
- 1.5308196 X 10-^°]?*meV, 



(55) 



where the constant term is again in perfect agreement 
with Borie's value 0.0050 meV [70] (Table 7). 



C. Finite size correction to the Kalian and Sabry 
contribution 



We can then evaluate the Kallen and Sabry contribu- 
tion AEjpj^^ - AE^^^^ using Vzi calciilated following Sec. 
Ill A, with good accuracy, using our numerical wave- 
functions. For a point nucleus I obtain 



ae: 



21,pn 

2si/2-2pi/2 



= 1.508097 meV, 



(56) 



in agreement with the resiilt of Ref. [37], 1.5079 meV, 
and in excellent agreement with the one from Ref. [39], 
1.5081 meV. Using the wavef unctions calculated with 



the proton size, I can also evaluate the finite size correc- 
tion to the Kallen and Sabry contribution. A direct fit to 
the numerical data gives a resiilt of the form 



AE^^'*" , (R) = 1.508097 - 0.00021341293]?^ 

2si/2-2pi/2^ ' 

+ 7.3404895 x W'^R^ (57) 
- 5.0291143 X 10-^]?*meV. 



and 



AEf , (]?) = 0.0000414300 - 9.25489 x IQ-^^R^ 

2pi/2-2p3/2 ^ ' 

+ 2.33622 X 10-"]?^ 
- 1.80063 xlO-i^K*meV. 



(58) 



for the fine structure, to be compared with 0.00004 meV 
in Ref. [70]. 



V. HIGHER-ORDER VACUUM POLARIZATION 
CORRECTIONS 

A. Higher-order vacuum polarization 

The term named "VP iteration", which correspond to 
Fig. 7, is given by Eq. (215) of Ref. [38] 

AEvPVp(2s) = 0.01244- - (Za)VrC^ (59) 

9 \7I/ 

where /jj. =94.96446 MeV for muonic hydrogen (using 
[11]). This adds 0.15086 meV to the Lamb-shift for 
muonic hydrogen. The Uehling potential under the form 
used in Sec.II can be introduced in the potential of the 
Dirac equation (3) when solving it numerically. This 
amounts to get the exact solution with any number of 
vacuum polarization insertions as shown in Fig. 7. The 
numerical methods that we used are described in Ref. 
[56, 57]. Because of the Logarithmic dependence of the 
point nucleus Uehling potential at the origin, we do not 
calculate the iterated vacuum polarization directly for 
point nucleus. We instead calculate for different mean 
square radii and charge distribution models, and fit the 
curves with /(]?) = a + bR^ + cR^ + rf]?*. All 4 models 
provides very similar values. The final value is 



^gii,ioop,fs ^ 0.15102161 

2Sl/2-2pi/2^ ' 



- 0.000098409804]?2 
+ 7.9038238 x IQ-^R^ 

- 7.2004764 X 10"^]?* meV. 



(60) 



The value calculated in Ref. [37] is 0.1509 meV and the 
one in Ref. [39] is 0.1510 meV in very good agreement 

with the present work. The method employed here pro- 
vides in addition the proton size dependence for this 



FIG. 7. Feynman diagrams obtained when the Uehling poten- 
tial is added to the nuclear potential in the Dirac equation. A 
double line represents a bound electron wavefunction or prop- 
agator and a wavy line a retarded photon propagator. The grey 
circles correspond to the interaction with the nucleus. Diagram 
(b) correspond to Fig. 5 and third term in Eq. (25) in Ref. [73] 



correction, which was not calculated before. For the fine 
structure, I obtain in the same way 



^^iijoop fs ^ 2.33197 X 10"^ 

- 1.77101 X IQ-'K^ 
9.36429 X 10-^°R3 
- 1.79951 xlO-^°K*meV. 



(61) 



FIG. 8. Lower order Feynman diagrams included in the Kal- 

len and Sabry Vxyij) potential, when the Uehling potential is 
included in the differential equation. See Figs. 4 and 7 for ex- 
planation of symbols. Diagrams (a) and (b) exactly correspond 
to diagrams (a) and (b) in Fig. 5 and Eq. (25) in Ref. [73] 



fine structure this correction is very small: 
^E2jS-j3/.(^) = 3-75754 xlO-« 



- 3.49318 X IQ-i^R^ 
-h 2.58244 X 10-"K^ 

- 2.74201 xlO-"R^meV. 



(63) 



B. Other higher-order Uehling correction 

Since we include the vacuum polarization in the Dirac 
equation potential, aU energies calculated by perturba- 
tion using the numerical wavefunction contains the con- 
tribution of higher-order diagrams where the external 
legs, which represent the wavefunction, can be replaced 
by a wavefunction and a bound propagator with one, 
or several vacuum polarization insertion. For example 
the Kallen and Sabry correction calculated in this way, 
contains correction of the t5^e presented in Fig. 8. This 
correction is given by 

A4'S-2;;aW = 0-0021552 



- 1.32976 X Vd'^'R} 

4- 9.4577 X 10-^r3 

- 8.5185 xlQ-^R^meV, 



(62) 



with a 10"^ meV accuracy. This correction is part of the 
three-loop corrections form Ref. [73, 74]. The diagrams 
in Fig. 8 correspond to diagrams (a) and (b) of Fig. 5 
in Ref. [73] and (e) (upper left) and (f) in Fig. 2 of 
Ref. [75]. The sum of contributions of the diagram (a) 
and (b) is 0.00223, in good agreement with our all-order 
fully relativistic result. The three loop diagram Fig. 5 
(c) Ref. [73] and Fig. 2 (g) Ref. [75] is included in the 
all-order contribution obtained by solving numerically 
the Dirac equation with the Uelhing potential. For the 



C. Wichmann and KroU correction 

We use the approximate potentials as presented in 
Refs. [68, 76] to evaluate the Wichmann and Kroll [77] 
Via correction to the Uehling potential. The correspond- 
ing diagram is shown on Fig. 4 (b). This contribution 
is given together with the light-by-light scattering dia- 
grams of Fig. 9 in Refs. [73, 75, 78]. We find for a point 
nucleus, the exact value, and a size correction, given by 

rl3,fs 



AE"''' , (R) = -0.0010170628 

2sl/2-2pi/2^ 



+ 5.5414179 X 10-^R2 

- 5.1356872 x lO'^^R^ 

- 1.9364450 X lO'^R* meV, 



(64) 



to be compared to the value given in Ref. [75] (Table III) 

of -0.001018(4) meV. In the lowest order approximation, 
the diagram on Fig. 9(a) provides an energy shift of 
AEis/Z^ [75, 78]. For the fine structure, this correction is 
comparable to the contribution from Eq. (63): 



^E2p5-2.3/.(^) = -4-21088 xlO-« 
4.41081 X 10-"r2 
- 8.49036 X 10-^5r3 
-M.81474X 10-^^2^* meV. 



(65) 
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FIG. 9. Feynman diagrams for the light-by-Ught scattering. See 
Figs. 4 and 7 for explanation of symbols. 




FIG. 10. Feynman diagrams for the muon self-energy. See Figs. 
4 and 7 for explanation of symbols. 



D. Muon radiative corrections 

1. Muon self-energy 

Highly accurate self-energy values for electronic 
atoms and point nucleus are known from Ref.[79]. The 
self-energy correction, represented in Fig. 10 is conve- 
niently expressed by the slowly varying function F(Za) 
defined by 



n ny 



(66) 



where m is the particle mass. 

The recoil corrections to F(Za) are described in detail 
in [11]. The dependence in the reduced mass has to be in- 
cluded leading to the following expressions, specialized 
for the n = l shells: 

a{Zaf ' ■■ 



AE 



fiSE,ZS 



n 8 

4 
+- 



(I) -1- 



4 10 
-lnfco{2S)+- 



na 



/67 

30 

In 



161n2\ 



3 / i a^^, 



a 



■G25(a)), 



(67) 



AE 



HSE,2pii2 



a {Zaf 



n 8 
103 
^180 ^\a^[i 



1 / OT, 



+ a^G2pv,(a) , (68) 



and 

AE 



a {Zaf 



n 8 
29 



m. 



+ - 



90 



ln(^-^^c^ + a^G2,,,Aa)^. (69) 



The Bethe logarithms are given by lnfco(2S) = 
2.811769893 and lnfco(2P) = -0.030016709 [80]. The 
the remainders are given by G2s(a) = -31.185150(90), 
Gzpipia) = -0.97350(20) and G2p,pia) = -0.48650(20) [79, 
81]. One then gets the exact muon self-energy for each 
state. For the 2s state, this gives 0.675150 meV instead of 
0.675389 meV. For the 2pi/2 I get -0.00916882 meV and 
for 2p3/2, 0.008393377 meV in place of 0.01424054 meV 
and -0.00332838 meV respectively, if one woiild use only 
the low order A40 term. 

The finite size correction is given by perturbation the- 
ory [11] Eq. (54) 

EsE-Ns{R,Za) = {A]n2 -jy(Za)SNs{R,Za), (70) 
where ([11] Eq. (51)) 

is the lowest-order finite nuclear size correction to the 
Coulomb energy. Here Ac = 1.867594282 fm is the muon 
Compton wavelength. Equation (71) provides Sns{R) - 
5.19745R^ for muonic hydrogen in agreement with Refs. 
[37, 39]. 

The self-energy correction to the Lamb shift with 
finite-size correction is then 



AE^J^^-AE^';'^'^'(R) = -0.68431882-^0.000824R^meV, (72) 
and to the fine structure: 



-SE,fs/ 



AE 



SE,f.s 
2p3/2 



AEf * = 0.017562197 meV. 



(73) 



It should be noted that in Ref. [39], the R^ -dependent 
part of the 2s self-energy is much larger than what is 
given in Eq. (72). This value was checked indepen- 
dently by using an all-order calculation with finite size, 
following the work of Mohr and Soff [82]. The results of 
this calculation agree reasonably well with Eq. (72) and 
is given by [83] 



-0.68431882 

-I-0.0012176389R2 
-0.00072582511R^ 
-I-0.00034744609R'' 
-0.000063788419R5meV. (74) 
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2. Muon loop vacuum polarization 

The vacuum polarization due to the creation of virtual 
muon pairs is represented by the same diagram 4 (a) 
and same equations (28) as vacuum polarization due to 
electron-positron pairs, replacing the electron Compton 
wavelength by the muon one. For S states, it is given by 
[11] Eq. (27),[35] Eq. (32) 



a{aZ)' 



\ 15 48/\m|, 



m^,(?, (75) 



in which higher order terms in Za have been neglected. 
For the 2s Lamb shift in muonic hydrogen it gives 
0.01669 meV and is included in Refs. [37, 39] and [84] 
Eq. (2.29) for the first a correction. As it is a sizable con- 
tribution, and the muon Compton wavelength, which 
represent the scale of QED corrections for muons is of 
the order of the finite nuclear size (1.9 fm), one could ex- 
pect a non-negligible finite size contribution. Using the 
numerical procedure described in Sec. 11, replacing the 
electron Compton wavelength by the muon one in Eq. 
(28), I obtain 

AEf (R) = 0.01671487464 

- 0.00005279721702R2 
+ 0.00001269866912R3 

(76) 

- 5.360546098 x lO'^K* 

+ 0.00001717649157R2 log(R) 

-H 2.047113814 X IQ'^R^ log(R) meV, 

where the constant term is in excellent agreement with 
(75) and the R dependence explicit. From Ref. [11], Eq. 
(55), one obtains 

^%p{ns) = ^a{Za)£Ns{R,a) = 0.000207581?^ meV (77) 

for the 2s level. This term is about 4 times larger than 
the numerical coefficient for in Eq. (76). 

Using the wavefunction evaluated with the Uehling 
potential in the Dirac equation, I also obtain the value of 
the sum of diagrams with one muon vacuum polariza- 
tion loop and any number of electron loops on each side, 
as in Fig. 7, with one loop being a muon loop: 

AE^^^!_l^'^^{R) = 0.00005348857 

- 5.358667885 x lO'^R^ 
+ 6.495888541 x lO'^K^ 



- 4.287455607 x IQ-^R* 
+ 2.822488184 x 10"^R^ log(R) 
+ 1.739448734 x lO-^RMog(R) 
meV. 



(78) 



AE: 



2pi/2-2p3/2^ 



This muonic vacuum polarization is a small contribution 
to the fine structure 

,,,iixii,fs^^^^ = 1.67794 x 10"^ 

- 2.10861 X 10-i°R2 
+ 8.51427 X 10-^^R^ 
- 1.38884 xlO-^iR^meV. 



(79) 



VI. EVALUATION OF THE RECOIL CORRECTIONS 

The relativistic treatment of recoil corrections is de- 
scribed in, e.g., [11], Eq. (10). The analytic solution of 
the Dirac equation for a point nucleus and a particle of 
mass m is given by 



Ed = mc^f(n, j) 



with 



f{n,i) = 



1 



1 + 



(80) 



(81) 



^ (n-;+i+V(i+|)'-(z«F) 
The recoil can then be included by evaluating [62, 85] 



Em = M(^ + iirC^{f{n,j)-l) 
2M 

l-5ijo (Za)V?c2 



(82) 



(83) 



k(2/ + 1) 2n3M2 

where M = m^i+ Mp. If one expands the previous equa- 
tion in power of (Za), one would find that the terms of or- 
der up to (Za)* are identical to what is given in Ref. [62]. 
We also compared the numerical results from our nu- 
merical approach for point nucleus, as described in Sec. 
I A to what can be obtained by using directly (80) and 
find excellent agreement. Below, we will make exclusive 
use of the direct numerical evaluation of the Dirac equa- 
tion. The relativistic corrections to Eq. (83) associated 
with motion of the nucleus are called relativistic-recoil 
correction. The correction to order (Za)^ and to all orders 
in m^/Mp is given by [11, 41, 85, 86] 

^ -^lnfco(U) 



"9" ~ 



7 



In 



X 



M^ln 



m. 



3 ""(Za)2 

Ml -ml 

2, Pp 
- In — 



(84) 



where 




/(/ + 1)(2/ + 1)' 



(85) 
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where 



FIG. 11. Feynman diagrams corresponding to the Relativistic 
recoil correction (84). The heavy double line represents the 
proton wave function or propagator. The other symbols are 
explained in Fig. 4. 



This correction corresponds to the diagrams in Fig. 11. 

The next order of the relativistic recoil corrections is 
given for s states by 



11 



In 



1 



6071 (Za)2 j' 
and for I > 1 states by 



(86) 



3- 



/(/ + 1) 



+ 3)}- 



(4/2-1) (2/ 
Using Eqs. (82) to (87), I obtain 



AE^f ,„ =0 + 0.0574706 - 0.0449705 
= 0.0125001meV, 
and to the fine structure: 

^^2py2-2p3,2 = 0.0000051 - 0.0862059 + 
= -0.0862008 meV. 



(87) 



(88) 



(89) 



This is in excellent agreement with the results from Ref. 
[70]. 



VII. EVALUATION OF SOME ALL-ORDER HYPERFINE 
STRUCTURE CORRECTIONS 

The expression of the h5^erfine magnetic dipole op- 
erator can be written as 



Hhfs = -ecoL ■ A{r) = -eca ■ A{r), 



with 



A(r) = 



/io X r 
An 



(90) 
(91) 



where /x is the nuclear magnetic moment and we have 
assumed a magnetic moment distribution of a point par- 
ticle for the nucleus. It is convenient to express Hiifg 
using vector spherical harmonics. On obtains [87-90] 



HHfs = M^-T\ 



(92) 



:a-r/;'(r) 



T\r) = 



-le- 



(93) 



and M.^ representing the magnetic moment operator 
from the nucleus. The operator acts only on the 
boimd particle coordinates. The vector spherical har- 
monic (r) is an eigenfunction of and /z, defined 
as [88, 90-93] 

Yf^ (f) = (r) = 2^ C (1, 1,1;^- o, o, q) Yi,,_, (r) 

a 

(94) 

where C(/i,/2//;'Wi,m2,m) is a Clebsh-Gordan coeffi- 
cient, Yi^ are scalar spherical harmonic and are eigen- 
vectors of and s^, the spin 1 matrices [88, 90-93]. The 
reduction to radial and angular integrals is presented 
in various works [88-90]. In heavy atoms, the hyper- 
fine structure correction due to the magnetic moment 
contribution is usually calculated for a finite charge dis- 
tribution, but a point magnetic dipole moment (see, e.g., 
[88, 89]). When matrix elements non-diagonal in / are 
needed, one can use [94] for a one-particle atom 



AE: 



HFS 
Ml 



2Mp Jo 



dr 



Pi{r)Q2{r) + P2{r)Qi{r) 



, (95) 



where g = jip/l = 2.792847356 for the proton, is the 
anomalous magnetic moment, A is an angular coefficient 



A = i-iy^i' 



+F 



I 

h 



111 

-I 1 



x(-l)^'4V(2/i + l)(2/2 + l)( 5 I i\]n{h,k,l2), 



(96) 



where n (h,k, h) = if /i + /i + 1 is odd and 1 otherwise. 
The ji are the total angular momentum of the i state for 
the bound particle, are orbital angular momentum, I 
is the nuclear spin, k the multipole order {k = 1 for the 
magnetic dipole contribution described in Eq. (95)) and f 
the total angular momentum of the atom. The difference 
between AE^^^ values calculated with a finite or point 
nuclear charge contribution is called the Breit-Rosenthal 
correction [95]. 

To consider a finite magnetic moment distribution, one 
uses the Bohr-Weisskopf correction [96]. The correction 
can be written [97] 



AE 



BW 



- -A^ r 

2Mp Jo 

r- , Pi 

X I dr— 
Jo 



drnrinirn) 
(r)Q2(r) + P2(r)Qi(r) 



(97) 
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where the magnetic moment density ju(r„) is normalized 

as 



f 

Jo 



dr„ri^{r„) = 1. 



(98) 



Borie and Rinker [98], write the total diagonal hj^er- 
fine energy correction for a muonic atom as 



AEy = 



47TK (F(F + 1)- 1(1 + l)-;0- + l)) ga 



X I dr 

JQ 



2M„ 



I dr„r^^BR{rn 
Jo 



)• (99) 



JO 

where the normalization is different: 

I d\^BRir„) = 471 I dr„rl^BR{rn) = 1- (100) 
Jo Jo 

This means that jUBR{r) = n{r)/{ATi). Evaluation of the 
Wigner 3J and 6J symbols in (96) give the same angiilar 
factor than in Eq. (95). 

The equivalence of the two formalism can be easily 
checked: starting from (99) and droping the angular fac- 
tors, we get, doing an integration by part 



r , Pi{r)Q2{r) r , 2 , 

J dr J dr„rf,^{r„ 

= ^J^ dr„ri^ir„) J dt— 



miit) 



Jo Jo 
Jo 



Pi{r)Q2{r) 



. fdry^Mr.) fdr^^, (101) 
Jo Jo 

where we have used (98). We thus find that the for- 
mula in Borie and Rinker represents the full hyperfine 
structure correction, including the Bohr-Weisskopf part. 

In 1956, Zemach [99] calculated the fine structure en- 
ergy of hydrogen, including recoil effects. He showed 
that in first order in the finite size, the HFS depends 
on the charge and magnetic distribution moments only 
through the Zemacs's form factor defined in Eq. (15). 
The proton is assumed to be at the origin of coordinates. 
Its charge and magnetic moment distribution are given 
in terms of charge distribution p(r) and magnetic mo- 
ment distributions jU(r). Zemach calculate the correction 
in first order to the hj^erfine energy of s-states of hy- 
drogen due to the electric charge distribution. The HFS 
energy is written as 



{Sp-s,) Jl 



(p{r) |2 i^{r)dr (102) 



(p the non-relativistic electron wavefunction and are 
the spin operators of the electron and proton. If the 



magnetic moment distribution is taken to be the one of a 
point charge, /./(r) - 6(r), the integral reduces to | (p{0) p. 
The first order correction to the wavefunction due to the 
nucleus finite charge distribution is given by 

cp{r) = (pc{0)ll-am^ J p(u)|tx - r|duj, (103) 

where (pciO) is the tmperturbed Coulomb wavefunction 
at the origin for a point nucleus. Replacing into Eq. (102) 
and keeping only first order terms, we get (Eq. 2.8 of 
Ref. [99] corrected for a misprint): 

X |l - lantfi p{u) \u-r \ (i{r)dvJr^ , 

= Ep |l - 2am^ p{u) \u-r \ ^{r)dudr^ , 

(104) 

where Ep is the well known HFS Fermi energy. Trans- 
forming Eq. (104) using r — > r + m, etc. Zemach obtains 



A£|ps = £F(l-2am^ (rz)). 



(105) 



with (rz) given in Eq. (16). The 2s state Fermi energy is 
given by 



^ 3 ''m„m 



(106) 



A. Hypetfine structure of the 2s level 

In order to check the dependence of the hyperfine 
structure on the Zemach radius and on the proton fi- 
nite size, I have performed a series of calculations for 
a dipolar distribution for both the charge and magnetic 
moment distribution. We can then study the dependence 
of the HFS beyond the first order corresponding to the 
Zemach correction. I calculated the h5rperfine energy 
splitting AEHFs(i^z,R) = Ehfs(R) + E^w (R,Rm) numer- 
ically. I also evaluate with and without self-consistent 
inclusion of the Uehling potential in the calculation, to 
obtain all-order Uehling contribution to the HFS energy. 
We calculated the correction AEhfs(Kz/R) for several 
value of Rz between 0.8 fm and 1.15 fm, and proton sizes 
ranging from 0.3 fm to 1.2 fm, by steps of 0.05 fm, which 
represents 285 values. The results show that the cor- 
rection to the HFS energy due to charge and magnetic 
moment distribution is not quite independent of R as 
one would expect from Eq. (105), in which the finite size 
contribution depends only on Rz. We fitted the h3^er- 
fine structure splitting of the 2s level, E?Lg (Rz, R) by a 
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hinction of R and Rz, which gives: 
E^pg(Rz,i^) = 22.807995 



- 0.0022324349K^ + 0.00072910794^^ 

- 0.000065912957^" - 0.16034434Rz 

- 0.00057179529Ri?z 

- 0.00069518048R2Rz 

- 0.00018463878R3Rz 
+ 0.0010566454i?2 
+ 0.00096830453RR2 

2r>2 



(107) 



+ 0.00037883473R^_R;^ 

- 0.00048210961R| 

- 0.00041573690RR^ 

+ 0.00018238754R^ meV. 

The constant term should be close to the sum of the Fermi 
energy 22.80541 meV and of the Breit term [100] . the HFS 
correction calculated with a point-nucleus Dirac wave- 
function for which I find 22.807995 meV. When setting 
the speed of light to infinity in the program I recover 
exactly the Fermi energy. The Breit contribution is thus 
0.002595 meV, to be compared to 0.0026 meV in Ref . [40] 
(Table II, line 3) and 0.00258 meV in Ref. [70]. Mar- 
tynenko [40] evaluates this correction, which he names 
"Proton structure corrections of order and a^", to be 
-0.1535 meV, following [35]. He finds the coefficient 
for the Zemach's radius to be -0.16018 meVf~^, in very 
good agreement with the present all-order calculation 
-0.16034 me Vf-i. Borie's value [70] -0.16037 meVf"! is 
even closer. The difference between Borie's value and 
Eq. (107) is represented in Fig. 12 as a function of the 
charge and Zemach radii. The maximum difference is 
around 1 jieV. 

In Ref. [35], the charge and magnetic moment dis- 
tributions are written down in the dipole form, which 
corresponds to (13), 



Gm [q^] 



A" 



l + Kp {A2+q2f' 



(108) 



with A = 848.5MeV. This leads to R = 0.806 fm as in 
Ref. [2] and Rz = 1.017 fm using this definition for 
the form factor in Eq. (17). Moreover there are re- 
coil corrections included. Pachucki [35] finds that the 
pure Zemach contribution (in the limit wZp oo) is 
-0.183 meV. In Ref. [101], the Zemach corrections is 
given as 6(Zemach) x Ep = -71.80 x lO'^Ef, for Rz = 
1.022 fm. This leads to a coefficient -0.1602 meVfm"^, in 
excellent agreement with our value -0.16036 meVfm"^. 

The effect of the vacuum polarization on the 2s hyper- 
fine structure energy shift as a function of the Zemach 
and charge radius have been calculated for the same set 
of values as the main contribution. The data can be de- 
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FIG. 12. difference between Borie's value and Eq. (107) result 
as a function of the charge and Zemach radii (meV). 



scribed as a function of Rz and R as 



Ehf7 (^Z'^) = 0.074369030 + 0.000074236132R2 

-I- 0.00013277334R3 - 8.0987285 x lO'^R'' 

- 0.0017880269Rz - 0.00017204505RRz 

- 0.00037499458R2Rz 

- 0.000070355379R3Rz 

- 0.00022093411R2 + 0.00035038656RR^ 
+ 0.00020554316R2R^ + 0.00025100642R3 

- 0.00017200435RR| 



0.000061266973R*meV. 



(109) 



It corresponds to the diagrams presented in Fig. 13. 
The size-independent term 0.07437 meV corresponds 
to the sum of the two contributions represented by 
the two top diagrams in Fig. 13 and is given as 
^^iioop-after-ioop VP = ^'^"^^^ ^^^V in Ref. [40]. The term 
AE^^yp = 0.0481 meV corresponds to a vacuum polar- 
ization loop in the HFS potential [35, 40, 98], which is not 
evaluated here. Corrections present in Ref. [70] not in- 
cluded in Eqs. (107) and (109) gives an extra contribution 
of 



eS£° = 0.10287meV. 



"HFS 



(110) 
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Combining Eqs. (107), (109) and (110), I get 

E^s (^Z'^) = 22.985234 - 0.0021581988R2 
+ 0.00086188128R3 

- 0.000074011685R* 

- 0.16213237Rz 

- 0.00074384033RRz 

- 0.0010701751R2]?z 

- 0.00025499415R3Rz (111) 
+ 0.00083571133R| 

+ 0.0013186911RR| 
+ 0.00058437789K^R| 

- 0.00023110319R| 

- 0.00058774125]IR| 

+ 0.00012112057K|meV. 

In Ref. [40], the equivalent expression is 

E^"^- (Rz) = 22.9857 - 0.16018Rz meV, (112) 
while it is 




FIG. 13. Feynman diagrams corresponding to the evaluation 
of the hyperfine structure using wavefunctions obtained with 
the Uehling potential in the Dirac equation. The grey squares 
correspond to the hyperfine interaction. 



equation: 



-■258006 

"HPS 



(Rz) = 22.9627 - 0.16037Rz meV, 



(113) 



in Ref. [70]. Using a Zemach's radius of 0.9477 fm in Eq. 
(112), needed to reproduce entry 11 in Table II of Ref. 
[40], one obtains 22.8148 meV as expected. In Eq. (113), 
it gives 22.8107 meV. Using the same Zemach radius 
and Eq. (Ill) I obtain 22.8104 meV with the muonic 
hydrogen proton radius value and 22.8103 meV with 
the COD ATA one, in excellent agreement with Borie's 
value. All three values are in agreement with the result 
22.8146(49) meV in Ref. [19]. In a recent work, how- 
ever, the use of Form factors in the Breit equations leads 
to smaller finite size corrections, leading to 22.8560 meV 
[102]. A comparison between some of these results is 
presented in Table I. 



VIII. EVALUATION OF MUONIC HYDROGEN n = 2 
TRANSITIONS 



(114) 



AE^°';Jl2pj^^(R) = 206.0209137 - 5.226135625R^ 
+ 0.03432100160R^ 
+ 0.0005454642475R'' 

- 0.00008785574420R^ 

+ 0.0002962967640R2 log(R) 

- 0.00004751147090R'' log(R) meV. 



In the same way Eqs. (46), (55), (58), (61), (63), (65), (73), 
(79), and (89) lead to the fine structure interval (which 
include the recoil corrections (89), included in Table II 
for the 2s Lamb shift) 



^£2;?-2,3/a(^) = 8-352051651 

- 0.00005203798087R2 

+ 1.215060759 x10-^R3 

- 1.544056441 x lO-'^R^'meV. 



(115) 



A. Lamb shift and fine structure 



The results presented in this work for the Lamb shift 

(Eqs. (45), (54), (57), (60), (62), (64), (72), (76), (78)) can 
be summarized in the following proton-size dependent 



Martynenko [105] finds E^°*'*'_2p^^^ = 8.352082 meV for 
the fine structure. 

A number of terms not included in Eqs. (114) are pre- 
sented in Table II together with the relevant references. 
Combining Eq. (114) with the sum of the contributions 
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TABLE I. Comparison of contributions to the 2s hyperfine structure from Refs. [40, 70] and the present work (meV) for 
Rz = 1.0668 fm [2] as used in Ref. [40]. Note that in this reference, the proton structure correction of order (item # 6) 
may combine the Zemach correction and the recoil correction (# 24). VP: Vacuum Polarization. 





# 


Ref. [40] 


Ref. [70] 


This work 


Fermi energy 


1 


22.8054 


22.8054 




Dirac Energy (includes Breit corr.) 


2 






22.807995 


Vacuum polarization corrections of orders , ct in 2nd-order 


3 


0.0746 


0.07443 




perturbation theory Cyvx 










All-order VP contribution to HPS, with finite magnetisation distribution 


4 






0.07244 


finite extent of magnetisation density correction to the above 

Proton structure corr. of order 


5 




-0.00114 




6 


-0.1518 


-0.17108 


-0.17173 


Proton structure corrections of order 


7 


-0.0017 






Electron vacuum polarization contributionH- proton structure corrections of order a* 


8 


-0.0026 






contribution of ly interaction of order 


9 


0.0003 


0.00037 


0.00037 


evplEp (neglected in Ref. [40]) 


10 




0.00056 


0.00056 


muon loop VP (part corresponding to £vP2 neglected in Ref. [40]) 


11 




0.00091 


0.00091 


Hadronic Vac. Pol. 


12 


0.0005 


0.0006 


0.0006 


Vertex (order a^) 


13 




-0.00311 


-0.00311 


Vertex (order a^) (only part with powers of ln(a) - see Ref. [103] ) 


14 




-0.00017 


-0.00017 


Breit 


15 


0.0026 


0.00258 




Muon anomalous magnetic moment correction of order a^, 


16 


0.0266 


0.02659 


0.02659 


Relativistic and radiative recoil corrections with 


17 


0.0018 






proton anomalous magnetic moment of order a* 










One-loop electron vacuum polarization contribution of ly interaction 


18 


0.0482 


0.04818 


0.04818 


of orders , a'' (eypi) 










finite extent of magnetisation density correction to the above 


19 




-0.00114 


-0.00114 


One-loop muon vacuum polarization contribution of ly interaction of order a* 


20 


0.0004 


0.00037 


0.00037 


Muon self energy +proton structure correction of order a* 


21 


0.001 




0.001 


Vertex corrections+proton structure corrections of order a* 


22 


-0.0018 




-0.0018 


"Jellyfish" diagram correction+ proton structure corrections of order a* 


23 


0.0005 




0.0005 


Recoil correction Ref. [104] 


24 




0.02123 


0.02123 


Proton polarizability contribution of order 


25 


0.0105 






Proton polarizability Ref. [104] 


26 




0.00801 


0.00801 


Weak interaction contribution 


27 


0.0003 


0.00027 


0.00027 


Total 




22.8148 


22.8129 


22.8111 



contained in Table II, I obtain the final 2s - 2pi/2 energy: 

= 206.0465137 - 5.2 

-I- 0.03432100160R3 



AE^°';f_2py^(R) = 206.0465137 - 5.226135625R^ 



(116) 



-I- 0.0005454642475]^* 

- 0.00008785574420]^^ 
+ 0.0002962967640]^2 ^^^^^^^ 

- 0.00004751147090]^'' log(]^) meV. 

This can be compared with the result from from E. Borie 
[70] 

^^Borieis ^ 206.0579(60) - 5.22713]?2 

2Sl/2-2pi/2^ ' V {117) 

0.0365(18)]?^ meV, 



and Carroll et al. [106] 
= 

-h 0.0546]?^ meV 



^ES/^-2pv2(^) = 206.0604 - 5.2794]?^ 



(118) 



An extra recoil contribution is given in Ref. [72] for 
the fine structure, corresponding to entry #9 in Table II 
for the Lamb shift. 



^4™-%3/. = -0.00006359 meV. 



(119) 



This term correspond to corrections beyond the fuU Dirac 
term. This lead to the final result 



^£S;f-2p3/.(^) = 8-351988061 



- 0.00005203798087]?2 

+ 1.215060759 x 10"^]?^ 

- 1.544056441 x IQ-'^R^meW. 



(120) 
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TABLE II. Contributions to the Lamb shift not included in Eq. (114) (meV). The uncertainty on the proton polarization value used 
in Ref. [12] has been increased by a factor of 10, according to the discussion in Ref. [26]. 



# 


Contribution 


Reference 


Value 


Unc. 


1 
2 

3 


NR three-loop electron VP (Eq. (11), (15), (18) and (23)) 

Virtual Delbriick scattering (2:2) 

Light by light electron loop contribution (3:1) 


[73] 
[75, 78] 

inn no^ 
[/D, /8J 


0.00529 
0.00115 
-0.00102 


0.00001 
U.UUOOl 


4 
5 


Mixed self-energy vacuum polarization 
Hadronic vacuum polarization 


[35, 84, 107] 
[108-110] 


-0.00254 
0.01121 


0.00044 


6 
7 

8 
9 


Recoil contribution Eqs. (82) and (83) 

Relativistic recoil of order (Za)^ Eq. (84) 

Relativistic Recoil of order (Za)* Eq. (86) 

Recoil correction to VP of order m/M and (m/M)^ in Eq. (4) 


[11, 36, 62, 85] 
[11, 37-39, 41] 
[11, 37] 
[72] 


0.05747063 
-0.04497053 
0.0002475 
-0.001987 




10 
11 


Proton Self-energy 
Proton polarization 


[35, 37,41,111] 
[18, 37, 109, 112, 113] 


-0.0108 
0.0129 


0.0010 
0.0040 


12 
13 


Electron loop in the radiative photon 

of order d'-l^af 

Mixed electron and muon loops 


[98, 114-116] 

[117] 


-0.00171 

0.00007 




14 


Rad. Recoil corr. a(Zdf 


[61] 


0.000136 




15 
16 

17 


Hadronic polarization a{Zafmr 
Hadronic polarization in the radiative 
photon a^XZafm-r 

Polarization operator induced correction 


[109, 110] 
[109, 110] 

[110] 


0.000047 
-0.000015 

0.00019 




18 


to nuclear polarizability a(Zafmr 
Radiative photon induced correction 
to nuclear polarizability a{Zafmr 


[110] 


-0.00001 






Total 




0.0256 


0.0041 



B. Transitions between hyperfine sublevels 



Using the results presented above I get 



ELt (^Z'K) - Eg;!, (Kz'K) = 209.92441 



The energies of the two transitions observed experi- 
mentally in muonic hydrogen are given by 



-f=2 

"2p3/2 



= AE: 



2si/2-2pi/2 + AE2pj/2-2p3/2 



, 3 p2p3/2 _ 1 p2s 

8 ™^ 4 HFS' 



(121) 



and 



^Iplll ^251/2 - ^E 



'2sl/2-2pi/2 + ^E2pi/2-2p3/2 

. 5 p2p3/2 ^ 3 25 f=l 
g^HFS 4^HFS ^ ''^HFS • 



(122) 



Here w^e use the results from [105] for the 2p states: 



Eots = 7.964364 meV 
E^f = 3.392588 meV 
6E^=P^ = 0. 14456 meV. 



(123) 



- 5.2261075^2 

+ 0.0343795061^3 
-I- 0.00043287446]^" 

- 0.000063788419R5 
+ 0.040533092Rz 

-1- 0.00018596008KKz 
+ 0.00026754376K2]?z 
-I- 0.000063748539K3Rz 

- 0.00020892783K| 

- 0.00032967277RRI 

- 0.00014609447i^22^| 
+ 0.000057775798R| 
+ 0.00014693531RR| 

- 0.000030280142R^ 

-I- 0.00029629676R2 \o^(R) 

- 0.000047511471R" log(R) 

(124) 

This can be compared with the result from 
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U. Jentschura [84] 



:.Jents.,F=2 



■ ^2^"^'^ = 209.9974(48) 
- 5.2262i?2 meV, 



(125) 



using the 2s hyperfine structure of Ref. [40]. 
For the other transition I obtain 

E^;; {Rz,R) - eL:° iRz,R) = 229.66162 



- 5.2282657R2 

+ 0.035241387R3 
+ 0.00035886278K^ 

- 0.000063788419]^^ 

- 0.12159928Kz 

- 0.00055788025KKz 

- 0.00080263129K2_Rz 

- 0.00019124562^3^2 
+ 0.00062678350K^ 

+ 0.00098901832i?K| 
+ 0.00043828342K2r^ 

- 0.00017332740K| 

- 0.00044080593KK^ 

+ 0.000090840426R| 

+ 0.00029629676R2 iog(|^) 

- 0.00004751 1471R* log(R) 

(126) 

Using Eq. (124), a Zemach radius of 1.0668 fm from 
Ref. [40] and the transition energy from Ref. [12], I 
obtain a charge radius for the proton of 0.84091(69) fm 
in place of 0.84184(69) fm in Ref. [12] and 0.8775(51) fm 
in the 2010 COD ATA fundamental constant adjustment. 
This is 7.1 a (using the combined a) from the 2010 CO- 
DATA value. A summary of proton size determinations 
is presented in Table III and Fig. 14. 



IX. CONCLUSION 

In the present work, I have evaluated finite-size de- 
pendent contributions to the n - 2 Lamb shift in muonic 
hydrogen, to the fine structure and to the 2s hyperfine 
splitting. The calculations were performed numerically, 
to all order in the finite size correction, in the framework 
of the Dirac equation. High-order size contributions to 
the Uelhing potential and to higher-order QED correc- 
tions been evaluated. The full dependance of the 2s hy- 



perfine splitting on the proton charge distribution and 
Zemach radius has been evaluated as well. 

The discrepancy between the proton size deduced 
from muonic hydrogen and the one coming from CO- 
DATA is slightly enlarged when tacking into account all 



0.91 



> 0.87 



0.81 
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FIG. 14. Plot of the proton size as a function of time and 
method. 



the newly calculated effects. It is changed from 6.9 a to 

7.1 CT. 



meV. 
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Appendix A: Coefficients for the numerical evaluation of 
the Kallen and Sabry potential for a point nucleus 

The functions defined in Eq. (35) are given here. We 
find, for x <3, the functions valid for a point nucleus: 

goir) = 0.00013575124407339550307r* 

- 0.00012633396034194731891r^ 
+ 0.00237541931 191 15541914r^ 

- 0.0052460271878852635132r^ 
4-0.16925588925254111005?-'' (Al) 

- 0.25201860708873574898?-3 
+ 0.95109984162919008905r^ 

- 2.0864972181198001792r 
-I- 1.6459704071917522632, 



gi{r) = -0.000078684672329473358699r^ 

- 0.0012293141869424835524r^ 

- 0.097906849416525020713r^ (A2) 

- 0.416662901899756662257^ 
-I- 0.13769050748433509769 

and 

g2(r) = -0.000012756169252850100497r^ 

+ 0.017425498169562658160r'' (A3) 
+ 0.44444444460943167625 . 

For X > 3, we fitted the coefficients in Eq. (40) to the 
ntraierical values. We obtain 

a = 4.3942926509010 
b = -10.059551479890 
c = 5.5493632222582 
d = 5.3327556570422 
e = -9.0762837836987 
/ = 5.1094977559523. 

Using these functions w^e reach an agreement to 9 deci- 
mal place with both the result of the numerical evalua- 
tion and the expansion from [68]. 
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Appendix B: Coefficients for the numerical evaluation of the and 
Kallen and Sabry potential for a finite nucleus. 

The coefficients for the functions defined in Eq. (39) 
that we obtained are listed below: 

ho(r) = -0.00001608988362060r' 
+ 0.000015791745043r* 

- 0.00036443366062r^ 
+ 0.0008743378646r^ 

- 0.038046259798r5 (Bl) 
+ 0.063004651772r^ 

- 0.363329158537-3 
+ 1.04324860906r2 

- 2.39716878893r + 2.005566300 

gi(r) = 0.7511983817345282548 
+ 0.13888763396658555408r2 
+ 0.020975409736870016795?-^ (B2) 
+ 0.00017561631242035479319r^ 
+ 9.057708511987165794 x lO'^r^ 



giir) = -OA 

- 0.0034850996339125316319r* (B3) 

- 1.4173521392055667219 x lO'^r* . 



